
% delete one (jack knife) estimate
% of impedance and variances
% Reference Chave, A. D., and Thomson, D. J., 1989
% J. Geophys. Res., 94, 14215-14225
% AG Jones (LiMS)
% Latest date 05.03.2003
% Minimizes the offdiagonal variances

% Tested.
%%%% 

function[tf9,tf_dof] = jknf(SPMatrix,ProcDef), % SPMatrix is a 4D array of cross & auto powers
% nfreq len nchanel nchannel

fprintf('jknf-> Jack knife estimate of tfs & var\n');
A = size(SPMatrix);
len=ProcDef.len;

for i = 1:A(1), % for all frequencies
    
   len=ProcDef.len;

   iflag = ones([len 1]);
   Z = jktf(SPMatrix(i,:,:,:),ProcDef,iflag);
   
   % June 24, 2011. To delete all un determined tfs
   
    L = any(isnan(Z) | isinf(Z),2);
    Z  = Z(~L,:);
    iflag(L) = [];
    len = len - sum(L);
   
   % end change 24 June 2011
   
   [Zjack,d1diff,var0] = jknfvar(Z,iflag);
   var1l = sqrt(var0(1)*var0(2)); % 2 -> Zxy 3-> Zyx
   j = 1;
   flag = 1;
   while j < (len-3) & flag == 1,
      kxy = find(d1diff(:,1) == max(d1diff(:,1)));
      kyx = find(d1diff(:,2) == max(d1diff(:,2)));
%       if length(kxy) > 1 | length(kyx) > 1,
%           disp('JKNF-> Warning length of kxy or kyx gt 1 using only first value');
%           disp(kxy);disp(kyx);
%       end;
      iflag([kxy(1) kyx(1)]) = 0; % kxy(1) = 0 as in some cases where length(kxy) > 1, there are problems 22.8.4
      n = sum(iflag);
      [Zjack,d1diff,var] = jknfvar(Z,iflag);
      var1 = sqrt(var(1)*var(2));
      if var1l > var1,
         var1l=var1;
         j = j+1;
      else
         iflag([kxy(1) kyx(1)]) = 1;% kxy(1) = 0 as in some cases where length(kxy) > 1, there are problems 22.8.4
         n = sum(iflag);
         [Zjack,d1diff,var] = jknfvar(Z,iflag);
         flag = 0;
      end;
   end;

tf9(i,1:4) = Zjack;
tf9(i,5:8) = var;
tf_dof(i) = n;  % d o f for this average
end;
ProcDef.len = n;